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AN ANALYTIC RECURSIVE METHOD FOR OPTIMAL MULTIPLE STOPPING: 
CANADIZATION AND PHASE-TYPE FITTING 


TIM LEUNG, KAZUTOSHI YAMAZAKI, AND HONGZHONG ZHANG 

Abstract. We study an optimal multiple stopping problem for call-type payoff driven by a spectrally 
negative Levy process. The stopping times are separated by constant refraction times, and the discount 
rate can be positive or negative. The computation involves a distribution of the Levy process at a constant 
horizon and hence the solutions in general cannot be attained analytically. Motivated by the maturity ran¬ 
domization (Canadization) technique by Carr [14], we approximate the refraction times by independent, 
identically distributed Erlang random variables. In addition, fitting random jumps to phase-type distribu¬ 
tions, our method involves repeated integrations with respect to the resolvent measure written in terms of 
the scale function of the underlying Levy process. We derive a recursive algorithm to compute the value 
function in closed form, and sequentially determine the optimal exercise thresholds. A series of numerical 
examples are provided to compare our analytic formula to results from Monte Carlo simulation. 
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1. Introduction 

A wide array of financial applications can be formulated as optimal multiple stopping problems. These 
include energy delivery contracts such as swing options [12, 13, 51], derivatives liquidation [25, 36, 37], 
real option analysis [15, 17, 19, 44], as well as employee stock options [24, 38, 39] potentially with 
additional reload and shout options [18]. In many of these applications, consecutive stopping times are 
separated by a constant or random period. In the literature, especially that of swing options, this timing 
constraint is commonly referred to as the refraction period. In real option analysis, the refraction period 
can be interpreted as the time required to build an infrastructure after an investment decision is made. 

In this paper, we discuss an analytic recursive method to solve a refracted optimal multiple stopping 
problem driven by a Levy process. This paper focuses on a computational aspect of the optimal multiple 
stopping problem. It is well known from related studies that the optimal strategy is of threshold-type. 
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Consequently, the optimal stopping problem reduees to finding these thresholds. However, the deter¬ 
mination of the threshold values still involves eomputing expeetations of a funetional at the end of the 
eonstant refraetion period, whieh is generally not explieit. In existing literature, Monte Carlo simulation 
methods are typieally employed to evaluate these expeetations (see [8, 45], among others). However, in 
praetiee this approaeh ean be eomputationally expensive and even infeasible in terms of the run time. 
Furthermore, with multiple stopping, one needs to know the entire expeeted future payoff funetional 
(with respeet to the starting point of the underlying proeess) in order to determine baekwards these fune- 
tionals as well as the optimal threshold levels for earlier stages. The simulation approaeh eommonly 
involves eomputing these expeetations for arbitrarily large number of starting points, and this adds to 
the eomputational burden and limits its applieability. In this regard, it is important to approximate these 
funetions in closed form so as to earry out effieiently the baekward induetion. 

One key feature of our analysis is that the rate for diseounting future eash flows ean be negative 
or positive. A negative diseount rate ean aeeommodate a number of applieations, sueh as stoek loans 
[11, 49] as well as real option problems where the investment eost grows faster than the risk-free rate. 
In these eases, one ean interpret that the effeetive diseount rate is negative. As argued by Blaek [9] (see 
also referenees therein), it is eommonly assumed that the nominal short rate must stay positive, but the 
real interest rate ean potentially be negative, espeeially during low-yield regimes. Henee, our framework 
permits diseounting eash flows at a negative effeetive or real interest rate. 

In our model, the underlying proeess is a speetrally negative Levy proeess, whieh has reeently been 
widely used in mathematieal finanee. Negative jumps ean model sudden downward movements of an 
asset priee. These proeesses are suitable in the struetural models of eredit risk and generate non-zero 
limiting value of the eredit spread as the maturity goes to zero as studied in [20, 26, 35, 40, 47]. Some 
reeent applieations of speetrally negative Levy proeesses inelude the prieing of perpetual Ameriean and 
exotie options [1, 6], optimal dividend problems [7, 33, 43], and eapital reinforeement timing [21]. For 
related optimal multiple stopping problems under speetrally negative models, we mention [51] for a 
swing put option with eonstant refraetion times, and [50] with a more general payoff funetion without 
refraetion times. For models with more general proeesses, Leung et al. [41] study a refraeted optimal 
multiple stopping problem driven by a two-sided Levy proeess with general random refraetion times, 
and Christensen and Lempa [16] eonsider a similar problem driven by a general Markov proeess with 
exponential refraetion times. 

Motivated by the maturity randomization (Canadization) method proposed by Carr [14], we provide 
an analytieal approximation by replaeing every eonstant refraetion time with an independent Erlang 
random variable, or a finite sum of independent, identieally distributed exponentially distributed times. 
Our method involves repeated integrations with respeet to the resolvent measure, whieh is written in 
terms of the scale function of the underlying speetrally negative Levy proeess. For the randomization 
methods applied in the prieing of finite-time horizon Ameriean options, we refer the reader to [28, 34]; 
similar ideas are also used in reeent work on the so-ealled Wiener-Hopf simulation [31, 23]. Bouehard 
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et al. [10] analyze a maturity randomization algorithm and apply it to stochastic control problems with 
applications to optimal single stopping and dynamic hedging under uncertain volatility. 

In order to apply the randomization method, the closed form expression must be preserved after the 
integration is applied with respect to the resolvent measure. This is satisfied when the Laplace exponent 
of the underlying Levy process has a rational form, in which case the scale function can be written as a 
finite sum of (possibly complex) exponentials (see [30]). Here, we focus on phase-type Levy processes 
[2], which constitute an important class of Levy processes with Laplace exponents of rational transform. 
In principle, any spectrally negative Levy process can be approximated by a Levy process of this form 
(which we call phase-type fitting). In particular, Egami & Yamazaki [22] give a series of numerical 
experiments for approximating the scale function of a general spectrally negative Levy process by that 
of a phase-type Levy process. For a hyperexponential fitting method applied to a CGMY process with a 
completely monotone Levy density, see [3]. 

Motivated by these, we combine phase-type fitting and randomization methods to compute efficiently 
the solutions of the optimal multiple stopping problem. Specifically, given a general spectrally negative 
Levy process, we first approximate it by a phase-type Levy process, and then approximate the solutions 
by randomizing the constant refraction times using independent, identically distributed Erlang random 
variables. We shall show that the resulting approximating value functions are written in closed form, 
with the associated parameters computed recursively. 

Our objective is to evaluate numerically the effectiveness of our approach, especially the accuracy of 
the value functions as a result of (i) phase-type fitting of the jump distribution, and (ii) refraction times 
randomization. Regarding part (i), while it is theoretically known that the class of phase-type distri¬ 
butions is dense in the class of all positive-valued distributions, there does not currently exist a single 
algorithm that can produce a sequence of phase-type distributions that are guaranteed to converge to a 
desired distribution (unless it has a completely monotone density). As for refraction times randomiza¬ 
tion, we refer to [10] and [42] for the related convergence results on the randomization approach. In a 
related study [28], detailed numerical experiments are conducted to confirm the convergence for pricing 
American put options when the Eevy process is in the meromorphic class. It is noted, however, that 
these results do not apply directly to our case, because we deal with a multiple optimal stopping problem 
where the refraction time is randomized. In addition, the payoff function (call type) is not bounded. For 
these reasons, it is important that our approach is evaluated numerically. 

In this paper, through a series of numerical examples, we show that our method is capable of accurately 
and efficiently computing the sequence of value functions and optimal exercise thresholds. In addition, 
the run-time analysis shows that this approach is significantly faster than the Monte Carlo simulation 
methods that adopt the Euler’s method to approximate the expected value function of the next stage at 
the constant refraction time. On the other hand, as the number of stages and the shape parameter of the 
Erlang distribution increase, the usual machine double precision may not be capable of computing the 
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parameters in the value funetions. Barring this potential issue, our elosed form formulas are eonfirmed 
by eomparing with the simulated values and allow for more effieient eomputation. 

The rest of the paper is organized as follows. Seetion 2 reviews the optimal multiple stopping problem 
for a speetrally negative Levy proeess and the eharaeterization of the optimal strategies in terms of up- 
erossing times. Seetion 3 presents our randomization method and derives the analytie value funetions 
reeursively using the resolvent measure. Seetion 4 diseusses the applieations of phase-type fitting and 
shows a baekward induetion formula to eompute the parameters of the value funetions for the randomized 
problem. We eonelude in Seetion 5 with numerieal evaluation of our proposed method. 

2. Preliminaries 

Let (fl, P) be a probability spaee hosting a speetrally negative Levy proeess X = {Xt)t>o- We 
define F := {Xt)t>o as the eompleted filtration generated by X, and T the set of all [0, cxo]-valued F- 
stopping times. We denote P^, as the probability and Ej. as the expeetation with initial value Xq = x G M. 
In partieular, when Xq = 0, we drop the subseripts in P^. and E^.. 

By the Levy-Khintehine formula, X ean be eharaeterized by its Laplace exponent given by 

(2.1) V'(s) := logE = cs + - 1 - S 2 rl{_i<^<o})n(d 2 r), s > 0, 

^ J (— 00 , 0 ) 

where If is a Levy measure with the support (—oo, 0) that satisfies the integrability eondition o) 
z^)T\.{dz) < oo. It has paths of bounded variation if and only if a = 0 and | 2 ;| n(d 2 ) < cxo; in this 
ease, we write (2.1) as 

(2.2) 'i/^(s) = cs+ [ (e^^^ — l)n(d 2 ), s > 0, 

J (— 00 , 0 ) 

with c := c — .2 n(dz). We exelude the ease in whieh —X is a subordinator (i.e., X has mono- 

tonieally deereasing paths a.s.). This assumption implies that c > 0 when X is of bounded variation. In 
addition, we assume throughout the paper that Xt admits a density; this is guaranteed to be satisfied if 
cr > 0 or the absolutely eontinuous part of the Levy measure has an infinite mass (see e.g. [48]). 

Assumption 2.1. We assume thatF{Xt G dx} dx for all t > 0. 

We eonsider the problem with sequential stopping (exereise) opportunities, where the payoff from 
eaeh exereise is 

(2.3) (t){x) :=e^ -K, xe M, 

for some eonstant K > 0. The assoeiated optimal multiple stopping problem is defined as 

' N 

,n=l 


v^^\x) := sup Ej- 


(2.4) 


X G M, iV G N. 
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Here, the optimization is over all inereasing sequenee of stopping times in sueh a way that any two eon- 
seeutive stopping times are separated by a refraetion period 5 > 0. In other words, the set of admissible 
strategies is given by 

(2.5) := {f = (tat, ... ,ri) e : r„+i + (5 < r„,n = iV - 1,..., 1}. 

Here we label in sueh a way that is the stopping time when there are n stopping opportunities left. 
Next, for any given diseount rate a G M, we define the proeess 

(2.6) :=Xt-at, f > 0, 


whieh is either a speetrally negative Levy proeess or the negative of a subordinator. As is well known, 
the limit of the running supremum X:= supo<t<oo a P-exponential random variable with rate 

parameter 

(2.7) $(a) := sup {A > 0 : '0(A) — Q!A<0} , 


with the eonvention that = oo a.s. when <h(a) = 0, and that X^^ = 0 a.s. when $(«) = oo. Henee, 
if <h(a) > 0, then the expeetation 

(2.8) ] = fl - ^<cx>, 

for any eonstant g e (0, $(«)). Sinee A i-A 0(A) — oA is strietly eonvex on [0, oo) and is zero at the 
origin, <h(a) > 1 if 0(1) < a. In this ease, we ean ehoose ^ > 1 sueh that the above moment generating 
funetion is finite, and thus, with the positivity of K, we have 


(2.9) 





< OO, X G M. 


This is a eritieal eondition so that the solution of the problem is nontrivial (see, e.g., [13, 41, 51]). In 
faet, we ean slightly weaken the eondition for the ease a < 0 (where exp{—at)K grows to infinity) to 
aeeommodate the ease 0(1) = a given 0'(1) < 0 (see [41] for a proof). 


Assumption 2.2. We assume that either (i) 0(1) < a or (n) 0(1) = a < 0 and 0'(1) < 0 holds. 


This guarantees that the value funetion is finite and admits a nontrivial solution. The optimal strategy 
is given by a sequenee of up-crossing times of the form 


( 2 . 10 ) 


* ._ 
Tn ■— 


r; := 1 o 0 .,. , +5 -f + 5, l<n<X-l, 


:= T+ o 

n+1 

for some parameters O* - io.n)l<n<N where 6 is the time-shift operator and 
(2.11) := inf {f > 0 : X* > a} , a G M, 
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with the usual convention that inf 0 = oo. The optimality of the threshold strategy for the single stopping 
problem has been shown by Mordecki [46]. The same characterization holds for the multi-stage problem, 
and we refer the reader to [13, 51] and the authors’ companion paper [41] for the proof. 

In view of these characterizations, the implementation of the optimal strategy reduces to identifying the 
values of a*, and this is the primary objective of our paper. To this end, we first rewrite (2.4) recursively 
as follows (see [29, 41], among others): 

(2.12) ^;(")(x) := supE4e-“"0(")(X.)l{,<^}] , 

where 

(2.13) (x) := ] , n = 1, 2,..., iV, 

and := 0. 

Given that an optimal stopping time in (2.12) is of threshold type, the value of a* can be determined 
by maximizing the value function over candidate threshold values: 

(2.14) a* G argmaxn^"'^(a:), 

agE 


(which maximizes uniformly in x G M), where 


(2.15) 


v^^\x) := 


'-{T+<oo} 


a, X G M. 


The following lemma is well-known for positive discount rate a > 0 (see [32], Theorem 3.12). Under 
Assumption 2.2, we generalize the result to accommodate the case with a < 0. Let 


(2.16) <h(a) := sup {A > 1 : 'ip{X) = a} , 

which is guaranteed to exist by Assumption 2.2 (which postulates that < a) and because xjj is 
strictly convex on [1, oo). 

Lemma 2.1. Under Assumption 2.2, we have l{r+<oo}] ~ 6xp(—$(a)|/)/or |/ > 0 and equals 

1 otherwise. 


Proof. We shall show for '0(1) < a (and hence $(«) > 1 by the strict convexity of -0); the case -0(1) = a 
then follows immediately by the monotone convergence theorem and the continuity of $(•). 

Fix 2 / > 0. The process (exp(<h(a)A'i — at))t>o is a martingale (see (3.11) of [32]), and hence we can 
derive (as in the first part of the proof of [32], Theorem 3.12) that 


(2.17) 


E 


<S>{a)X^^^+-a{tAT+) 

6 y 


1, t > 0. 


Here the integrand of the left-hand side is bounded in t by an integrable random variable, i.e.. 


(2.18) 
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where the last inequality holds beeause <l>(a) > 1 and X^^rp+ < y a.s. due to the lack of positive jumps. 
Hence applying dominated convergence in (2.17) gives 


(2.19) 


1 = E 




4t+<oo} 




e-«r+ ^ 


{T+ <oo} 


where the last equality holds as Xrp+ = y on {T+ < oo}. This completes the proof. 


□ 


Due to Lemma 2.1 and to the fact that the process X necessarily creeps upward and hence Xj.+ = a 
on {T+ < oo} under Pj. with x < a, we can write 


( 2 . 20 ) 




e-‘*’(“)(“-^)0(^)(a), x<a, 
X > a. 


Remark 2.1. It can be shown that the threshold levels are bounded from below by \ogK and increase 
as the number of remaining stopping opportunities decreases, i.e., log K < < • • • < a}. It has 

been shown in [41] that this monotonicity also holds when the refraction times 6’s are generalized to be 
independent, identically distributed random variables provided that they are independent of X, and Xs 
admits a density. They also show that there exists a limit a*^ := lim 7 v->.oo Ojv — log 


3. Recursive Analytic Formula 

The characterization of the optimal strategy as described in the previous section greatly simplifies the 
problem. In practice, however, the solution cannot be obtained analytically because in general the distri¬ 
bution of X^ is not known in view of (2.12). The biggest hurdle therefore is to compute the expectation 

(3.1) E^ [e-^^v^^-^\Xs)] , 2 <n<N. 

In order to circumvent this difficulty, we adopt the Canadization technique by Carr [14] and approximate 

(3.1) by replacing the constant 6 with some independent Erlang random variable r]{M, A), or equivalently 
a sum of M independent, identically distributed exponential random variables with parameter A. Herein, 
we set 

(3.2) A = Al^l := M/5. 

Then, r]{M, A) ~ 5 for large M by the strong law of large numbers. In other words, we solve the 
randomized version of the optimal multiple stopping problem in order to approximate the one with 
constant refraction times. For optimal stopping problems with random refraction times, the filtration 
needs to be modified; see [16] for the precise construction. However, this technical detail does not affect 
the resulting threshold structure of the optimal stopping strategies, as discussed in [16, 41]. 

In this section, we shall show that the value functions with randomized refraction times can be obtained 
recursively via the resolvent measure written in terms of the scale function. 
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3.1. First Step. We first construct the base case. In view of (2.20) for n = 1, because = 0, the 
value of al is obtained analytically via (2.14). The first order condition becomes 

(3.3) 0 = 4>'{a) — $(a)0(a), 


which admits a unique solution given by 
(3.4) al 


log 


$(a)/f 

$(«) - 1 ’ 


It is easy to check that this is equivalent to the smooth fit condition n(^)'(a^+) = ). 

We now start at 


(3.5) 




4>{x), X > 

X < a*, 


and derive an analytical expression for the expectation 


(3.6) 




as an approximation of (3.1) for n = 2. 

The very initial task is to compute the case of exponential time horizon, 

poo 

(3.7) := E, = A / dt = 

Jo 

where we define, for any measurable /, that (whenever it exists) 

(3.8) MJ:=X [ 0("+“)(x,d|/)/(|/) 

Jr 

for a resolvent measure 

poo 

(3.9) 0('')(a:, d?/) := / {X* G dy} df, y eR and g > 0. 

Jo 

It is known for the case of spectrally negative Levy process that this resolvent measure admits a density 
and can be written in terms of the so-called scale function. Fix g > 0, the (g-)scale function, 

(3.10) M ^ [0,CX)), 


is zero on (— 00 , 0), continuous and strictly increasing on [0, 00 ), and is characterized by the Laplace 
transform: 


(3.11) 


e-^^W^'i){x)dx = 


f{s) - q 


s > $(g), 


where 


( 3 . 12 ) 


<I)(g) := sup{A > 0 : '0(A) = g}. 
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By Corollary 8.9 of [32], the resolvent measure 0 ^'^) {x, dy) has a density (y — x) with respeet to the 
Lebesgue measure where 

(3.13) := ^ e M. 

Henee on eondition that A + a > 0, (3.7) ean be rewritten using the seale funetion. 

This manipulation can be applied repeatedly by further adding more independent, identically dis¬ 
tributed exponential random time horizons. Indeed, for any 2 < m < M and x G M, 

poo 

(3.14) m(1’”^)(x) = A / 

Jo 


3.2. Multiple Steps. Now that (3.1) is approximated by u^^’^\x) we can approximate (p^'^\x) as in 
(2.13) by 

(3.15) (j)^^\x) := (j){x) +u^^’^\x). 


Using this approximation, we can obtain an approximation to a^, say a^, by the first order condition 
0 = 0*^^^'(a) — <I>(a)0*^^^(a), and obtain an approximation to 


(3.16) 


U^‘^’^\x) = v^‘^\x) 


(j)^'^\x), X > a 2 , 

0(2)(a*)e-^0A52-0, a;<a^. 


Similarly to the first step, for any 1 < m < M, 


(3.17) u^^’^\x) := E, 




5 


which gives an approximation for (3.1) for n = 3. 

Continuing in this fashion, we can derive the approximations defined by 


(3.18) 

(3.19) 

(3.20) 

(3.21) 


4 >^'^\x) := 4 >{x) + ^’^\x), 

a* G arg{a G M : 0^"'^ (a) — <I)(a)0*'”^(a) = 0}, 

0(”)(a:), ^ x>a*^, 

0(”)(a;)e-‘^(“ASn-0^ x<al, 

uM^x) := E, 1 < m < M. 


u 


in, 0 )^x) =v^^\x) : = 


Finally, (x) is the desired approximation to our multiple stopping problem. For the rest of the paper, 
we let a] := a] for notational convenience. 


10 


T. LEUNG, K. YAMAZAKI, AND H. ZHANG 


4. Spectrally Negative Phase-type Case 


In order to carry out the algorithm described in the previous section, it is important that the backward 
induction can be done analytically. That is to say, the closed-form expression must be preserved after the 
operator AA as in (3.8) is applied. In this section, we shall show that this is possible if we focus on the 
phase-type Levy process of the form (4.1) below. 

It is known from Proposition 1 of [2] that, for any spectrally negative Levy process X, there exists a 
sequence of spectrally negative phase-type Levy processes converging to X in the space oo) 
of real-valued right-continuous functions with left-limits (cadlag); this implies that Xi in distri¬ 

bution (see Remark 1 of [2] and Corollary VII 3.6 of [27]). From this, we naturally conjecture that the 
resolvent of X can be approximated by that of a phase-type Levy process. Indeed, Egami & Yamazaki 
[22] show numerically that the scale function of a general spectrally negative Levy process can be ap¬ 
proximated at least when the Levy measure is finite. For the case of infinite Levy measure satisfying 
the condition of Asmussen and Rosinski [5], Asmussen et al. [3] show that a Brownian motion can be 
used as a proxy to approximate the frequent infinitesimal jumps (where they consider a Levy measure 
with a completely monotone density). In the next section, we analyze the approximation errors through 
numerical experiments. 

Throughout this section, let X be a spectrally negative phase-type Levy process. 



(4.1) 


0 < f < oo, 


for some c G M and a > 0. Here B = {Bt)t>o is a standard Brownian motion, N = {Nt)t>o is a Poisson 
process with arrival rate p, and Z = {Zn)n=i, 2 ,... is an independent, identically distributed sequence of 
phase-type-distributed random variables with representation {d, a, T). These processes are assumed to 
be mutually independent. Recall that a distribution on (0, cxo) is of phase-type if it is the distribution of 
the absorption time in a finite state continuous-time Markov chain consisting of one absorbing state and 
d G N transient states. Thus, any phase-type distribution can be represented by d, the d x d transition 
intensity matrix over all transient states T, and the initial distribution of the Markov chain a. 

Let t be the transition probabilities from the d transient states to the absorbing state. The Laplace 
exponent (2.1) is then 


(4.2) 


■0(5) = cs -f -I- p {cy.{sl — T) — l) , 


which can be extended to s G C except at the negative of eigenvalues of T. 
For the rest, let us define 


(4.3) 


p := a + X = a + M/5, 


and we assume this to be strictly positive. 
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Suppose {—i G Xp} is the set of the roots of the equality -0(5) = P with negative real parts. As 
is discussed in Section 5.4 of [30] and has been confirmed numerically in [22] (see also our numerical 
results in Section 5 of this current paper), it is highly unlikely that any root in Xp has multiplicity larger 
than one. Hence, we can assume that the roots in Xp are distinct. Consequently, the scale function can be 
written 


(4.4) 


W^P\x) 


0, a; < 0, 


where 

(4.5) 


• — 


5 + 


p-i}{s) 




see [22]. Here i G Xp} and [Ki^p] i G Xp} are possibly complex-valued. Hence the resolvent density 
(3.13) is written 


(4.6) 


e^p\z) 


$'(p)e--J-(p)^^ > 0, 

EiGX, ^<0. 


Our objective here is to show that the function for each n > 1 and 0 < m < M (derived 

recursively as in the previous section) is a piecewise function with subdomains (fo*, (see 

Remark 2.1 regarding the monotonicity of a*) where we define, for notational convenience, Oq := oo 
and := —oo. More specifically, we shall show that, on each subdomain, it is a sum of products of 
polynomials and exponentials: 


(4.7) 




for a* < X < where we define 

^n,m 


(4.8) 


+ ^ 1 < / < n + 1, 1/ G M, 


h=0 


with m := (n — l)M + m — 1. Here the parameter set 

(4.9) 

:= GXp,0 < h < < h < 

satisfies 

^(n,m,n+l) _ ^{n,m^n-\-l) _ ^(n,m,n+l) _ g 


(4.10) 
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Its proof and the derivation of the parameter set ^ can be done inductively. Along the same line 
as the arguments in the last section, we go through the backward induction. First, the base case (n = 1 
and m = 0) is trivial because, in view of (3.5), the function can be written as (4.7) by setting 

yl(i,o,i) = ^( 1 , 0 , 1 ) _ ^( 1 , 0 , 2 ) _ exp(-$(a)a*) with 4,^ = -1. 

In view of our discussion in the previous section, there are two types of inductive steps. The first 
kind increments the step counter n while the second kind increments m by applying the integration with 
respect to the resolvent measure. We shall call the former Step I and the latter Step II. 


4.1. Inductive Step I. We show that if the hypothesis holds for n > 1 and m = M, then it also holds 
for some n + 1 and m = 0. By this hypothesis, the equations (3.18) and (4.7) give, for a* < x < ai_^. 


(4.11) -K) + + l)e^+ 




i&Xp h=0 


h=0 


and 


(4.12) 0(”+i)'(x) = + l)e^ + Ci, 


h ' 




3 


X 


+ he 


^h-l\ 


X 


iGT-p h=0 


^n,M 

h=0 


By (3.19), we can identify the optimal threshold 

Now, in view of (3.20), the representation of ^(’^+^’0) can be obtained by setting 

(4.13) ^(n+l,0,n+2) _ 


and yl(«+i.o,n+2) ^ ^{n+i,o.n+2) ^ fjin+i,o,n+2) ^ ^(n+i,o,n+2) = Q, and for 1 < / < n + 1 


(4.14) 


^(n+l,0,0 _ ^ 


and 


Q{n+l,0,l) _ ^{n,M,l) jj{n+l,0,l) _ 

It can be confirmed that because, by assumption, have = 

^(n+1,0,1) ^ g_ 
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4.2. Inductive Step II. It is now sufficient to show that, for fixed n > 1 and 0 < m < M — 1, the 
hypothesis for m implies that for m + 1. 

For every x G M, with 1 < L(= L^) < n + 1 be sueh that a\ < x < a\_i. 


(4.16) 

T 


= 1 


{L>2} L 

KKL-l 




J X 

'^°'Li^Xp L+l<l<n+l''°‘’l i£Xp 

+ ^ (al, X, -ii^p) + l{L<n} ^ 


i&Ip 


where 

(4.17) 


:= I e ™/^”’”"’‘^(i/)di/, s<t 


ieip 


—qy r{n,m,l) / 


L+l<l<n+l 


In partieular, for x > a* (or L = 1), 


(4.18) 

A 


= oo, <h(p)) 


i&Xp 


i&Ip 


2<l<n+l 


while, for x < a* (or L = n + 1), 


= e*"’>-1>'(p) i(p)) + e*''’>">I>'(p)03i”+7>(i,<, <»(?)) 


(4.19) 


A 


(n,m) I 


l<l<n 

{n,m)i 




i&Xp 


By repeatedly applying integration by parts, the right-hand side of (4.16) ean be written in a elosed 
form in the form (4.7). As the eomputation is straightforward but tedious, we defer the proof of the 
following proposition and the detailed expressions of the reeursive formula (4.20) and the parameter set 
to the appendix. 
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Proposition 4.1. Fix n > 1 and 0 < m < M — 1 and assume (4.7), (4.9) and (4.10) hold. Then, with 
L = Lx being the unique integer such that dl<x< dl_ We have 


(4.20) 




A 


A 


4n,m4-l 




/7 = 0 

for some parameter set 
(4.21) 

Tn,m+1 ■■= i eXp,h> 0}, > Q|^ ^(n,m+l,0j 

that satisfies (4.10). 


l<l<n+l 


5. Numerical Results 


In this section, we evaluate our method numerically. For X, we use a spectrally negative Levy process 
of the form (4.1) with a modification so that (Z„) is independent, identically distributed sequence of the 
following random variables: 

Case 1: Exponential random variable with parameter 1, 

Case 2: Weibull random variable with parameter (2,1), 

Case 3: The absolute value of the Gaussian (folded normal) with mean zero and variance 1, 
whose respective densities are 


(5.1) 


exp{ —x}, 2x exp { —} 




X G (0, cxo). 


Case 1 is a special case of a phase-type random variable and its scale function can be computed exactly 
in the form (4.4). For Cases 2 and 3, the EM-algorithm is applied to approximate them by phase-type 
distributions for d = 6: the fitted phase-type distributions are 


’ -5.6546 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 ’ 


' 0.0000 ' 

0.6066 

-5.6847 

0.0000 

0.0166 

0.0089 

5.0526 


0.0007 

0.2156 

4.3616 

-5.6485 

0.9162 

0.1424 

0.0126 


0.9961 

5.6247 

0.0000 

0.0000 

-5.6786 

0.0000 

0.0000 

, (y. 

0.0000 

0.0107 

0.0000 

0.0000 

5.7247 

-5.7420 

0.0000 


0.0001 

0.0136 

0.0000 

0.0000 

0.0024 

5.7022 

-5.7183 


0.0031 
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and 


' -4.0488 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 ' 


' 0.0052 ' 

0.1320 

-4.0012 

0.0000 

0.0455 

3.7040 

0.0044 


0.0659 

0.2367 

0.8595 

-4.2831 

0.1897 

0.2918 

2.3724 


0.7446 

3.1532 

0.0000 

0.0000 

-4.0229 

0.0000 

0.0000 

, ol 

0.0398 

0.2497 

0.0000 

0.0000 

3.7024 

-4.0124 

0.0000 


0.0043 

0.0434 

2.1947 

0.0938 

0.1704 

0.1217 

-4.9612 


0.1403 


for Cases 2 and 3, respectively. For this phase-type fitting, we use EMpht which is written in C and is 
publicly available\ For more details of this method, we refer to [4]. Here, we choose d = 6 because, as 
has been confirmed in [22], the fitting can be conducted accurately and quickly. While accuracy tends to 
increase in d, run time increases nonlinearly; hence d cannot be chosen arbitrarily large. As we shall see 
below, our choice of d attains very small fitting errors. 

For numerical illustration, we set the negative discount rate a = —0.02 along with K = 100. For each 
case of Z, we consider the Levy process with common parameters p = 1.5 and a = 0.2 and choose 
so that '0(1) = a — 7 (i.e. exp(—(a — 7 )t + Xt)t>o is a martingale) for our choice of 7 . Notice that 
(and hence EXi as well) decreases as 7 increases. In the context of stock loans, as discussed in [11, 49], 
the negative discount rate is the difference of the risk-free rate and the loan rate, K is the loan amount, 
and 7 is the dividend rate. All the numerical results given below are generated by MATLAB scripts with 
double precision on a Windows 7 computer with an Intel Xeon CPU E5 — 2620, 2.00GHz, 24.0GB RAM. 


5.1. One-stage randomization. We first analyze the accuracy and computation time of our randomiza¬ 
tion algorithm by considering the expectation, for 5 = 0.5, 

(5.4) 

where is analytically given as in (3.5). In order to do so, we evaluate the approximations by our algo¬ 
rithm in comparison to the simulated results. More specifically, we first compute, for M = 1,..., 5,10, 
the approximations by these two methods to 

(5.5) 

and then approximate the constant 5 case (5.4) by simulation with a starting point x = a\. This enables 
us to analyze the approximation errors of our randomization algorithm for the Erlang case (5.5), and also 
analyze how large M needs to be to acquire accurate approximations for the constant 5 case (5.4). 

Eollowing the arguments in the previous section, our computation involves two main steps: (i) com¬ 
puting the roots of 0(-) = p, and (ii) computing recursively the parameter set F in (4.9). The root-finding 
procedure is conducted by MATLAB built-in function solve (). In Table 1, we give sample values 

'Available at http;//home.imf.au.dk/asmus/pspapers.html as of March 14, 2014. 
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7 = 

0.02 





7 = 

0.1 




M = 

1 

M = 

3 

M 

= 

1 

M = 

3 

,p 

1.0252 

- 1 - 

O.OOOOi 

1.5941 

- 1 - 

O.OOOOi 

1.0056 

- 1 - 

O.OOOOi 

1.5825 

- 1 - 

O.OOOOi 

^2,p 

3.8602 

- 1 - 

3.6058i 

3.9134 

- 1 - 

3.3255i 

3.8296 

- 1 - 

3.6319i 

3.8939 

- 1 - 

3.3384i 

^3,p 

3.8602 

- 

3.6058i 

3.9134 

- 

3.3255i 

3.8296 

- 

3.6319i 

3.8939 

- 

3.3384i 

^4:,p 

7.8211 

- 1 - 

3.4389i 

7.6518 

- 1 - 

3.2454i 

7.8398 

- 1 - 

3.4933i 

7.6613 

- 1 - 

3.2799i 

^5,p 

7.8211 

- 

3.4389i 

7.6518 

- 

3.2454i 

7.8398 

- 

3.4933i 

7.6613 

- 

3.2799i 

^6,p 

9.5837 

- 1 - 

O.OOOOi 

9.3632 

- 1 - 

O.OOOOi 

9.6386 

- 1 - 

O.OOOOi 

9.3983 

- 1 - 

O.OOOOi 

^7,p 

42.040 

- 1 - 

O.OOOOi 

46.026 

- 1 - 

O.OOOOi 

38.4292 

- 1 - 

O.OOOOi 

42.666 

- 1 - 

O.OOOOi 


Case2; Weibull 





7 = 

0.02 





7 = 

0.1 




M = 

1 

M = 

3 

M = 

1 

M = 

3 

^i,p 

0.9842 

- 1 - 

O.OOOOi 

1.4669 

- 1 - 

O.OOOOi 

0.9674 

- 1 - 

O.OOOOi 

1.4583 

- 1 - 

O.OOOOi 

^2,p 

3.2497 

- 1 - 

2.3023i 

3.2876 

- 1 - 

2.0887i 

3.2331 

- 1 - 

2.3200i 

3.2784 

- 1 - 

2.0976i 

CO 

3.2497 

- 

2.3023i 

3.2876 

- 

2.0887i 

3.2331 

- 

2.3200i 

3.2784 

- 

2.0976i 

'b4,p 

5.5298 

- 1 - 

1.62971 

5.4233 

- 1 - 

1.54371 

5.5425 

- 1 - 

1.6464i 

5.4300 

- 1 - 

1.5543i 

^5,p 

5.5298 

- 

1.6297i 

5.4233 

- 

1.54371 

5.5425 

- 

1.6464i 

5.4300 

- 

1.55431 

^e,P 

6.4520 

- 1 - 

O.OOOOi 

6.2947 

- 1 - 

O.OOOOi 

6.4805 

- 1 - 

O.OOOOi 

6.3103 

- 1 - 

O.OOOOi 

^7,p 

37.565 

- 1 - 

O.OOOOi 

41.862 

- 1 - 

O.OOOOi 

34.049 

- 1 - 

O.OOOOi 

38.617 

- 1 - 

O.OOOOi 


Case 3: Folded Normal 

Table 1. Values of for M = 1,3 and 7 = 0.02,0.1 (listed in ascending order). 
We can confirm that these values are all distinct. Because X has a Brownian motion 


component, |Xp| = d + 1 = 7. 


of here we can confirm that these values are all distinct. The step (ii) can be done efficiently by 
applying, for M times. Inductive Step II in the previous section. As for the simulated results, we com¬ 
pute this via Monte Carlo simulation based on 1 million sample paths, where the Brownian motions are 
approximated by random walks with time step At = T /100 for each inter-arrival time T between jumps. 

Tables 2 and 3 summarize the results for 7 = 0.02 and 0.1, respectively. The functions ^(i,M)>s, as 
in (5.5), obtained from the analytic recursive formula and simulation are listed for M = 1,..., 5,10, 
along with the constant 5 case computed by simulation presented in the bottom row. We also report 
the computation times (in seconds). The times that correspond to the analytic formula are given as a 
sum of the time spent for steps (i) and (ii). For the values under simulation, we give the mean and 95% 
confidence interval for each case. 

The simulated results are subject to some errors arising from the discretization of Brownian motions, 
but they are useful as a benchmark. These discretization errors are confirmed to be minimal in view of 
the comparison between these two methods for Case 1. Recall that the numerical results for Case 1 by 
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randomization 

simulation 

M 

value 

time 

value 

time 

1 

1823.65 

0.306 + 0.008 

1823.89(1821.61,1826.17) 

146.712 

2 

1824.27 

0.242 + 0.024 

1824.15(1822.03,1826.28) 

152.969 

3 

1824.51 

0.245 + 0.066 

1824.58(1822.51,1826.66) 

157.001 

4 

1824.64 

0.240 + 0.146 

1823.69(1821.71,1825.68) 

162.082 

5 

1824.72 

0.238 + 0.271 

1825.10(1823.00,1827.19) 

167.068 

10 

1824.88 

0.273 + 2.008 

1823.11(1821.01,1825.20) 

186.893 

const 

N/A 

1823.90(1821.80,1826.00) 

141.833 


Case 1: Exponential 



randomization 

simulation 

M 

value 

time 

value 

time 

1 

1665.62 

0.604 + 0.036 

1665.68(1663.64,1667.73) 

395.063 

2 

1666.12 

0.455 + 0.201 

1663.54(1661.63,1665.44) 

393.945 

3 

1666.32 

0.374 + 0.599 

1666.02(1664.07,1667.97) 

404.590 

4 

1666.42 

0.524 + 1.350 

1664.73(1662.84,1666.62) 

403.772 

5 

1666.49 

0.534 + 2.520 

1665.47(1663.59,1667.34) 

411.322 

10 

1666.58 

0.403 + 18.74 

1667.07(1665.08,1669.06) 

426.856 

const 

N/A 

1666.61(1664.75,1668.47) 

386.804 


Case 2: Weibull 



randomization 

simulation 

M 

value 

time 

value 

time 

1 

1482.88 

0.594 + 0.036 

1486.05(1484.40,1487.69) 

141.351 

2 

1483.35 

0.380 + 0.232 

1484.31(1482.76,1485.86) 

147.091 

3 

1483.53 

0.389 + 0.584 

1484.30(1482.72,1485.89) 

152.59 

4 

1483.63 

0.379 + 1.302 

1484.06(1482.45,1485.67) 

156.25 

5 

1483.69 

0.382 + 2.471 

1485.29(1483.69,1486.89) 

161.923 

10 

1483.80 

0.329 + 18.64 

1485.24(1483.67,1486.81) 

184.049 

const 

N/A 

1485.35(1483.78,1486.91) 

137.375 


Case 3: Folded Normal 

Tab LE 2 . Comparison between results under randomization and simulation for 7 = 0.02. 
The comparison is done for each Erlang shape parameter M ; in the bottom row (labeled 
const), the approximated values under simulation for the constant 5 = 0.5 case are given. 
The listed values under simulation are the mean and 95% confidence interval. The com¬ 
putation times (in seconds) for randomization are given as a sum of the time spent for 
steps (i) and (ii). 
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randomization 

simulation 

M 

value 

time 

value 

time 

1 

323.83 

0.360 + 0.008 

323.85(323.44,324.27) 

146.511 

2 

324.33 

0.235 + 0.023 

324.10(323.69,324.51) 

150.866 

3 

324.54 

0.231 + 0.064 

324.23(323.82,324.64) 

157.942 

4 

324.65 

0.235 + 0.143 

324.13(323.74,324.53) 

162.298 

5 

324.72 

0.237 + 0.273 

324.51(324.10,324.91) 

168.973 

10 

324.87 

0.271 + 2.013 

324.48(324.06,324.90) 

184.726 

eonst 

N/A 

324.97(324.56,325.37) 

142.774 


Case 1: Exponential 



randomization 

simulation 

M 

value 

time 

value 

time 

1 

303.13 

0.307 + 0.036 

302.36(301.94,302.78) 

389.083 

2 

303.54 

0.271 + 0.202 

303.44(303.04,303.85) 

396.285 

3 

303.72 

0.588 + 0.679 

303.63(303.23,304.03) 

401.650 

4 

303.81 

0.377+ 1.458 

303.48(303.09,303.86) 

402.040 

5 

303.87 

0.401 + 2.504 

303.65(303.24,304.06) 

404.430 

10 

304.00 

0.335 + 18.69 

303.87(303.48,304.25) 

428.611 

eonst 

N/A 

303.98(303.60,304.37) 

385.145 


Case 2: Weibull 



randomization 

simulation 

M 

value 

time 

value 

time 

1 

265.46 

0.876 + 0.037 

265.67(265.33,266.00) 

144.601 

2 

265.85 

0.380 + 0.235 

265.90(265.57,266.23) 

150.759 

3 

266.01 

0.386 + 0.583 

266.37(266.04,266.70) 

155.104 

4 

266.10 

0.341 + 1.311 

266.49(266.17,266.81) 

158.389 

5 

266.15 

0.388 + 2.459 

266.69(266.37,267.01) 

162.124 

10 

266.28 

0.336 + 18.75 

266.50(266.18,266.81) 

182.452 

eonst 

N/A 

266.86(266.55,267.17) 

139.149 


Case 3: Folded Normal 

Table 3. Comparison between results under randomization and simulation for 7 = 0.1. 
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the randomization algorithm are exact in the sense that there is no approximation error from fitting the 
seale function. Based on this observation, we can also infer from the results on Cases 2 and 3 that the 
associated fitting errors of the scale function are also minimal. This suggests the practicability of the use 
of the phase-type distribution as an approximation for a general Levy process. 

As M increases from 1 to 10, the approximate value function (5.5) increases monotonically and ap¬ 
proaches the simulated value for (5.4) associated with the constant 5 case. In fact, the exponential 
refraction time case (i.e. M = 1 ) already gives a reasonable approximation. 

In terms of the computation time, the randomization method is significantly faster than simulation. 
Note also that, for the randomization method, this computation is required only once to obtain the whole 
shape of the value function. The simulation method, on the other hand, is unfortunately not practical; 
it takes several minutes to attain this accuracy for a particular point of x. Recall that in our multiple 
stopping problem, we need to know the whole shape to conduct backward induction. If the simulation 
method is applied, one needs to compute for arbitrarily large number of starting points x. However, this 
is computationally infeasible. 

While the randomization method runs instantaneously when M is small, we observe that the compu¬ 
tation time increases nonlinearly in M. It also depends on the number of phases; Case 1 (with 1 phase) 
runs faster than Cases 2 and 3 (with 6 phases). This suggests one limitation of the randomization algo¬ 
rithm that the value of M and the number of phases d cannot be chosen arbitrarily large. However, as we 
already see in Tables 2 and 3, the approximate value function stabilizes even for small M and our choice 
of d. 

5.2. Multiple-stage case. We now move on to the multiple-stage case. Using our randomization algo¬ 
rithm, the approximate value functions , ■ ■ ■, are computed for Erlang shape parameters M = 1, 3 
and are shown in Figures 1 and 2, respectively, for 7 = 0.02 and 0.1. The threshold levels dl,... ,al 
(circles) are marked on the approximate value function curves. In particular, the top curve corresponds 
to the approximate value function As expected, the thresholds are all above the strike iT = 100 and 
they admit the ordering <+1 < a*. This is consistent with Remark 2.1. 

Recall that the process exp(—(a — 7 )^ + Xt)t>o is a martingale under the given parameters. Hence, 
for a small value of 7, the value function is close to linear in exp(x). On the other hand, as 7 increases, it 
appears to be more convex. Moreover, the function decreases as 7 increases because 7 reduces the 
drift of X. As in the single stopping case, the difference between the value functions for M = 1 and 3 is 
close to invisible. This suggests that these are reasonable approximations for the constant 6 case. On the 
other hand, the optimal threshold levels show non-negligible difference between the cases M = 1 and 3. 

5.3. Dependence on N and M. In Figure 3, we show the threshold levels with respect to the number 
of stages N and to the Erlang shape parameter M based on Case 3 with 7 = 0.1. On the left panel, 
we plot a^,..., 04 for fixed M = 1,..., 4. Note that the first threshold dl is independent of M. In the 
example on the right panel, we plot the threshold over M = 1,..., 10. As M increases from 1 to 10, 
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exp(x) 



exp<x) 


Case 1 (Exponential) with M = 1 


Case 1 (Exponential) with M = 3 



exp(x) 



exp<x) 


Case 2 (Weibull) with M = 1 


Case 2 (Weibull) with M = 3 




exp(x) exp<x) 

Case 3 (Eolded Normal) with M = 1 Case 3 (Eolded Normal) with M = 3 


Eigure 1. The approximate value funetions when 7 = 0.02 with threshold levels 
al,... ,al (eireles) marked on the approximate value funetion eurves. The values are 
monotone in the number of stages (the top eurve eorresponds to the approximate value 
funetion 
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Case 1 (Exponential) with M = 1 


Case 1 (Exponential) with M = 3 




Case 2 (Weibull) with M = 1 


Case 2 (Weibull) with M = 3 




Case 3 (Eolded Normal) with M = 1 


Case 3 (Eolded Normal) with M = 3 


Eigure 2. The approximate value funetions when 7 = 0.1. 
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the threshold first deereases relatively fast within the narrow range (5.81, 5.82) for M = 1, 2, 3, and then 
flattens toward the value 5.805 for larger M. Between M = 9 and 10, the difference is well less than 
0 . 001 . 

Figure 4 illustrates for iV = 1,..., 15 with 7 = 0.05 and M = 1. With more remaining exercise 
opportunities (large N), the function is higher and the optimal threshold for the previous exercise is 
lower. We observe that the distance between successive optimal thresholds (marked by circles) reduces 
as the number of remaining exercises increases (see e.g. the top value function curve marked with 15 
circles). 




Figure 3. Dependence of the thresholds on N and M for Case 3 with 7 = 0.1. The left 
panel plots a*,..., 04 for fixed M = 1, ..., 4. The right panel plots the threshold over 
M = 1,...,10. 

5.4. Limitations. Recall from formula (4.7) for recovering the function from the parameter set 
r = {A, B,C, D, E). In particular the coefficients D’s are multiplied by exp(<I)(a + M/5)x)x^, so 
this term tends to become very large near a*, even though H’s are zero above al- From our numerical 
tests, it can take value up to the order of 10^° while the values of D’s tend to remain small. Recall that 
p := a + M/5 increases in M and so does $(p) = <I)(a + M/5). In addition the maximum value of 
h (the counting index in (4.7)) increases as M and N increase. As a result, the computation can break 
down when the Erlang shape parameter M and/or the number of exercises N are large. MATLAB or 
other softwares with double precision cannot handle the computation involving these large numbers. 

In Figure 5, we plot the function computed by MATLAB for Case 3 with 7 = 0.1 and = 5 
for M = 4 (left) and M = 5 (right). While the parameters in F can be computed instantly and do not 
explode, numerical imprecision in computing the value function may arise when small parameters are 
multiplied by very large numbers and summed up. Indeed, with M = 4, discontinuities appear between 
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exp(x) 


Figure 4. Value functions and optimal thresholds (marked by circles) for Case 3 with 
= 1,..., 15 when M = 1 and 7 = 0.05. 


X 10 




exptx) exp(x) 

Figure 5. Limitations: the value functions computed by MATLAB for Case 3 with 
7 = 0.1 for M =4 (left), 5 (right). 

0,2 and al, and with M = 5 the error becomes visibly clear, yielding an inaccurate value function and 
threshold levels a*. This is consistent with the observation given in [28] (that deal with an American 
put option), where their randomization algorithm requires more than double precision. This issue can 
potentially be resolved by setting the machine epsilon as in [28] so as to increase precision. However, 
this is beyond the scope of our paper, because it requires special skills in computer science and it is our 
aim to evaluate the performance that can be achieved in a usual computing environment. Even given this 
potential limitation, the analytic formula is useful in its own right as it reveals the mathematical structure 
of the solution to the optimal multiple stopping problem. 
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This observation also highlights the potential trade-off between seleeting large values of M and N 
given limits on machine precision. Nevertheless, we have seen from Tables 2 and 3 that the approxima¬ 
tion remains stable for different small values of M. In summary, the analytic formula (4.7) is very useful 
and tractable for solving the optimal multiple stopping problem, as compared to the simulation approach. 


Appendix A. Proof of Proposition 4.1 and the updating formula 

Fix n > 1 and 0 < m < M — 1 and suppose (4.7), (4.9) and (4.10) hold. We shall show that the 
identity (4.16) can be written as (4.20) where the parameter set T^ ^+1 is given by (A. 10) below. 

First, for a < 13, straightforward integration gives the following expression for as in (4.17). 
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(2) When q = equals 






^i,p 






ii,p + 1 


E 

9=0 


- S®+^) 


9+^ 


Mi,p-ij,p)^ oh 'sp ^(n,m,0_ 9^ _ 

j&Tp\{i} h=Q g=h 




tL ^ t t 0 t 0^ t It 

(A.3) “ X] ^ Qirhrn,l) 


j&Ip\{i} h=0 

In.m 


g=h 


h\{Q,p-^^,p)P+^ 


—h 


Mi,p+'^(p))^ 
j=0 h=j 

ete,p+‘i’(“))i _ gte,p+^(“))« 


h\ 


L(«.,P+^(p))^oi _ p{^i,P+Hp))tfj \ _ 


-j 


_|_ £{n,m,l) _ 


^i,p + *^*( 0 ) 


In particular, if s = —00 and = Q{n,m,i) _ (j{n,m,i) _ 


(A.4) 


In,m I'i , 

t, -e.,p) = - 


h\ 


j=0 


h=j 




_|_ ]^in,m,l) 


e(Ci.P+'^(«))I 
^i,p + *^*( 0 ) 


By letting s = x and t = and multiplying by exp($(p)a;), we obtain 


(A.5) e'^^P^^WL{x,al_^, <I>(p)) = + B^n,m,L)^. 


In,m +1 


in.m +1 


E E + E (A”’”'- 


h)^^(p)x^h-^ ^{n,m,L) ^<S>{a)x ^ 


iGTp h=0 


h=0 


where 


^{n,m,L) ._ Q{n,m,L) ._ Biri,m,L)! 


(A.6) 


/^{n,m,L) 
'^i,h 


:=Eci 


{n,m,L) 

9 


0 < h < In,my 


ClrLY^O, i€lm 


g=h 














26 


T. LEUNG, K. YAMAZAKI, AND H. ZHANG 


and 
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with, for all i G X 
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Substituting (A.5) and (A.8) in (4.16), Proposition 4.1 is satisfied by setting X,m-i-i = X,m + 1 and 
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